/*
Copyright (c) [2019年5月1日] [吴超]
[MBT_studio] is licensed under Mulan PSL v2.
You can use this software according to the terms and conditions of the Mulan PSL v2.
You may obtain a copy of Mulan PSL v2 at:
		 http://license.coscl.org.cn/MulanPSL2
THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
See the Mulan PSL v2 for more details.
*/
#include "PBF_Simulation.h"

#include "core/CUDA_utils.cuh"


#define DEF_CUDA
#include <线性代数/几何计算.h>
#include <线性代数/矢量计算.h>
//#include <math.h>


//#define DEF_帧率 0.015
#define DEF_帧率 0.05
static float32 g防止除0补偿值 = 0.00000001f;


//初始化外力和加速度
__global__ void f_forceand(
	vec3* 位置,
	vec3* 上一帧坐标,
	vec3* 速度,
	vec3* 加速度,
	vec3* 外力,
	float32*	重量,
	//float32*	平方范数,
	//S_DomProp	域属性,
	vec3		重力 

	)
{
	//uint64 id = blockIdx.x + (blockIdx.x * gridDim.x) + (gridDim.x * gridDim.y * blockIdx.z);
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	
	上一帧坐标[当前ID] = 位置[当前ID];
	//printf(" 速度 %f %f %f =  %i \n", 外力[当前ID].x, 外力[当前ID].y, 外力[当前ID].z);
	if (重量[当前ID] != 0.0f) {
		//其它力施加
		加速度[当前ID] = vec_add(外力[当前ID], 重力);
		//加速度[当前ID] = 重力;
		//printf(" 速度 %f %f %f =  %i \n", 速度[当前ID].x, 速度[当前ID].y, 速度[当前ID].z);
		速度[当前ID] = vec_add(速度[当前ID], 加速度[当前ID] * DEF_帧率);
		位置[当前ID] = vec_add(位置[当前ID], 速度[当前ID] * DEF_帧率);
	}
	//printf(" 外力 %f %f %f =  %i \n", 外力[当前ID].x, 外力[当前ID].y, 外力[当前ID].z);
	vec_multi(&(外力[当前ID]), 0.95);

	//移动变化量[当前ID] = { 0,0,0 };
	//位置[当前ID] = vec_Max(位置[当前ID], vec_add(域属性.m_域世界偏移, 域属性.m_单元大小*2));
	//位置[当前ID] = vec_Min(位置[当前ID], vec_add(域属性.m_域世界偏移, vec_sub(域属性.m_世界大小, 域属性.m_单元大小*2)));
	//位置[i] = vec_Max(位置[i], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移);
	//位置[i] = vec_Min(位置[i], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移 + m_邻接.m_Data.m_域属性.m_单元属性.m_世界大小);
	//平方范数[id] = vec_平方范数(加速度[id] * 初始更新时间步长 + 速度[id]);
	//printf(" key x =  %f", vec_平方范数(vec_add(速度[id], vec_multi(加速度[id], 初始更新时间步长))));
	//printf(" key x =  %f", vec_平方范数(速度[id] + 加速度[id] * 初始更新时间步长));
	//平方范数[id] = vec_平方范数(速度[id] * 初始更新时间步长);
	return;
	//平方范数[id] = 0.512;
	//vec3 v = vec_add(速度[id], 加速度[id]) * 初始更新时间步长;
	//平方范数[id] = vec_平方范数(v);
}


//清空值 准备附加外力数据
__global__ void f_ClearValue(vec3* 移动变化量, vec3 重力) {
	//uint64 id = blockIdx.x + (blockIdx.x * gridDim.x) + (gridDim.x * gridDim.y * blockIdx.z);
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	移动变化量[当前ID] = 重力;

}


__global__ void f_Clampdown(vec3* 位置, S_DomProp 域属性) {
	//uint64 id = blockIdx.x + (blockIdx.x * gridDim.x) + (gridDim.x * gridDim.y * blockIdx.z);
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	位置[当前ID] = vec_Max(位置[当前ID], vec_add(域属性.m_域世界偏移, 域属性.m_单元大小 * 2));
	位置[当前ID] = vec_Min(位置[当前ID], vec_add(域属性.m_域世界偏移, vec_sub(域属性.m_世界大小, 域属性.m_单元大小 * 2)));
	
}




__device__ uvec3 f_getCoord_To_DomCoord(S_DomProp& dom, const vec3& coord) {
	return { uint32(coord.x / dom.m_单元大小), uint32(coord.y / dom.m_单元大小), uint32(coord.z / dom.m_单元大小) };
}


__device__ S_DomUnitIndex f_getUnitParticleIndex(S_cuAdjacencyIndex& adjacency, S_DomProp& dom, const uvec3& coord, vec3 c = {0,0,0}, uint64 _id = 0) {
	uvec3 分辨率;
	分辨率.x = dom.m_分辨率.x;
	分辨率.y = dom.m_分辨率.y;
	分辨率.z = dom.m_分辨率.z;


	uint64 id = f_坐标转线性索引(coord, 分辨率);
	//if (id >= 27000000) {
	//	printf("id =  %llu \n", id);
	//	printf("f_PBF_Gradient 位移步进 x =  %u %u %u \n", coord.x, coord.y, coord.z);
	//	printf("f_PBF_Gradient c x =  %f %f %f %llu \n", c.x, c.y, c.z, _id);
	//}
	S_DomUnitIndex Unit;
	//printf("adjacency.m_IndexNum =  %i \n", adjacency.m_IndexNum);
	//return {};
	Unit.num = adjacency.m_IndexNum[id];
	//printf("Unit.num =  %i \n", Unit.num);
	//return {};
	id = adjacency.m_VoxelIndexs[id];
	//printf("adjacency.m_VoxelIndexs = id %i \n", id);
	Unit.index = &(adjacency.m_Indexs[id]);

	return Unit;
}

__device__ vec3 f_PBF_Gradient(vec3 dir, S_cuAdjacencyIndex& 邻接) {
	vec3 梯度 = {};

	float32 rl = vec_normalize(&dir);
	const float32 q = rl / 邻接.m_粒子光滑半径;
	//printf(" rl =  %f", rl);
	//printf(" f_PBF_Gradient q =  %f", q);
	//printf(" f_PBF_Gradient rl =  %f", rl);
	//printf(" f_PBF_Gradient 邻接.m_L =  %f", 邻接.m_L);
	//if (q > 1.0f || rl < 0.00006f) {
	//	return {};
	//}
	if (q <= 1.0f && rl > 0.00006f) {
		
		vec3 q梯度 = vec_multi(dir, (1.0f / (rl * 邻接.m_粒子光滑半径)));
		if (q <= 0.5) {
			梯度 = vec_multi(q梯度, (邻接.m_L * q * (3.0f * q - 2.0f)));
		}
		else {
			const float32 系数 = 1.0 - q;
			梯度 = vec_multi(q梯度, (邻接.m_L * (-系数 * 系数)));
		}
		
	}
	//printf("f_PBF_Gradient 位移步进 x =  %f \n", 梯度.x);
	//printf("f_PBF_Gradient 位移步进 y =  %f \n", 梯度.y);
	//printf("f_PBF_Gradient 位移步进 z =  %f \n", 梯度.z);
	return 梯度;
}

__device__ float32 f_PBF_DensityWeights(vec3 dir, float32 m_K, float32 粒子光滑半径) {
	float32 rl = vec_normalize(&dir);
	const float32 q = rl / 粒子光滑半径;

	float32 权重 = 0.0f;
	if (q <= 1.0 && rl > 0.000001) {
		if (q <= 0.5) {
			const float32 q2 = q * q;
			const float32 q3 = q2 * q;
			权重 = m_K * (6.0f * q3 - 6.0f * q2 + 1.0f);
		}
		else {
			权重 = m_K * (2.0f * powf(1.0f - q, 3));
		}
	}

	return 权重;
}




__global__ void f_getUnitParticle(S_cuAdjacencyIndex adjacency, S_DomProp dom, uvec3 coord, vec3* pos, vec3* particles, uint32* num) {
	//uvec3 coord = f_浮点坐标转索引坐标(x, dom.m_域世界偏移, dom.m_单元大小);
	//__shared__ S_DomUnitIndex Unit;
	//if (threadIdx.x == 0) {
	//printf("f_getUnitParticle ======================= \n");
	S_DomUnitIndex Unit = f_getUnitParticleIndex(adjacency, dom, coord);
	(*num) = Unit.num;
	//}
	//__syncthreads();
	//printf(" key x =  %i", Unit.num);
	for (uint8 j = 0; j < Unit.num; ++j) {
		uint64 k = Unit.index[j];
		//printf(" Particle ID =  %i \n", k);
		particles[j] = pos[k];
	}
	//uint64 k = Unit.index[threadIdx.x];
}




__global__ void f_SolvingDensity(
	S_CP_parameter cCP,
	//vec3* 位置,
	float32* 质量,
	float32* c密度,

	float32* 拉格朗日,
	//uint32* c物理体ID,
	//uint8* c物理体类型,

	float32 密度0,
	float32 m_权重0,
	S_cuAdjacencyIndex	邻接索引,
	S_DomProp			域属性
)
{
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	vec3* 位置 = cCP.m_位置;

	vec3 x = 位置[当前ID];
	//float32 密度 = c密度[当前ID];
	float32 Density = 质量[当前ID] * m_权重0;
	//printf("m_权重0 =  %f \n", m_权重0);
	const static float32 eps = 1.0e-6;
	uint32 当前粒子物理ID = cCP.m_物理体ID[当前ID];
	uint32 当前粒子物理类型 = cCP.m_物理体类型[当前ID] & 0b1;

	uvec3 coord = f_浮点坐标转索引坐标(x, 域属性.m_域世界偏移, 域属性.m_单元大小);
	//return;
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 搜索坐标;
				搜索坐标.x = coord.x + ox;
				搜索坐标.y = coord.y + oy;
				搜索坐标.z = coord.z + oz;


				S_DomUnitIndex 单元索引 = f_getUnitParticleIndex(邻接索引, 域属性, 搜索坐标, x, 当前ID);
				//printf("单元索引.num =  %i \n", 单元索引.num);
				for (uint8 j = 0; j < 单元索引.num; ++j) {
					
					uint64 k = 单元索引.index[j];
					//uint32 邻接粒子物理ID = c物理体ID[k];
					//uint8 因子 = uint8(邻接粒子物理ID == 当前粒子物理ID);
					//if (k > 5000) {
						//printf(" 因子 ID =  %i \n", 因子);
					//}
					Density += f_PBF_DensityWeights(vec_sub(x, 位置[k]), 邻接索引.m_K, 邻接索引.m_粒子光滑半径);
					//printf("Density =  %f \n", Density);
				}
			}
		}
	}
	c密度[当前ID] = Density;
	//printf("Density =  %f \n", Density);


	const float32 约束 = DEF_Max(Density / 密度0 - 1.0, 0.0);
	float32 梯度和 = 0.0f;
	vec3 邻接梯度 = {};
	//const float32 约束 = DEF_Max(Density / 密度0, 0.0);
	//printf("Density =  %f \n", Density);
	//printf("约束 =  %f \n", 约束);
	if (约束 != 0) {
		for (int8 oz = -1; oz < 2; ++oz) {
			for (int8 oy = -1; oy < 2; ++oy) {
				for (int8 ox = -1; ox < 2; ++ox) {
					uvec3 搜索坐标;
					搜索坐标.x = coord.x + ox;
					搜索坐标.y = coord.y + oy;
					搜索坐标.z = coord.z + oz;


					S_DomUnitIndex 单元索引 = f_getUnitParticleIndex(邻接索引, 域属性, 搜索坐标);
					//printf("单元索引.num =  %i \n", 单元索引.num);
					for (uint8 j = 0; j < 单元索引.num; ++j) {
						uint64 k = 单元索引.index[j];
						//const float32 约束 = DEF_Max(密度 / 密度0 - 1.0, 0.0);
						uint32 邻接粒子物理ID = cCP.m_物理体ID[k];
						uint8 因子 = !(uint8(邻接粒子物理ID == 当前粒子物理ID) & 当前粒子物理类型);
						//printf(" 因子 ID =  %i \n", 因子);
						//(-质量[k] / 密度0);
						vec3 梯度 = f_PBF_Gradient(vec_sub(x, 位置[k]), 邻接索引);
						//printf("质量[k] =  %f \n", 质量[k]);
						梯度 = vec_multi(梯度, (-质量[k] / 密度0) * 因子);
						//printf(" 质量[k] / 密度0 =  %f \n", 质量[k] / 密度0);
						//printf(" 位移步进 x =  %f %f %f %i \n", 梯度.x, 梯度.y, 梯度.z, 因子);
						//printf(" 位移步进 x =  %f \n", 梯度.x);
						//printf(" 位移步进 y =  %f \n", 梯度.y);
						//printf(" 位移步进 z =  %f \n", 梯度.z);
						梯度和 += vec_平方范数(梯度);
						//printf("梯度和 =  %f \n", 梯度和);
						vec_sub(&邻接梯度, 梯度);
					}
				}
			}
		}


		梯度和 += vec_平方范数(邻接梯度);
		拉格朗日[当前ID] = -约束 / (梯度和 + eps);
	}
	else {
		拉格朗日[当前ID] = 0.0f;
	}
	
	//printf("约束 =  %f \n", 约束);
	//printf("梯度和 =  %f \n", 梯度和);
	//printf("拉格朗日[当前ID] =  %f \n", 拉格朗日[当前ID]);
}

	 
__global__ void f_ProjectionDelta(
	S_CP_parameter cCP,
	//vec3* 位置,
	vec3* 位置量变化,
	//float32* 密度,
	float32* 质量,
	float32* 拉格朗日,

	float32 密度0,
	S_cuAdjacencyIndex	邻接索引,
	S_DomProp			域属性

) {
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	uint32 当前粒子物理ID = cCP.m_物理体ID[当前ID];
	uint32 当前粒子物理类型 = cCP.m_物理体类型[当前ID] & 0b1;

	vec3* 位置 = cCP.m_位置;
	vec3 x = 位置[当前ID];

	//vec3 位移步进 = 位置量变化[当前ID];
	vec3 位移步进 = {};

	uvec3 coord = f_浮点坐标转索引坐标(x, 域属性.m_域世界偏移, 域属性.m_单元大小);
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 搜索坐标;
				搜索坐标.x = coord.x + ox;
				搜索坐标.y = coord.y + oy;
				搜索坐标.z = coord.z + oz;

				S_DomUnitIndex 单元索引 = f_getUnitParticleIndex(邻接索引, 域属性, 搜索坐标);
				//printf(" f_ProjectionDelta 单元索引.num =  %i \n", 单元索引.num);
				for (uint8 j = 0; j < 单元索引.num; ++j) {
					uint64 k = 单元索引.index[j];

					vec3 邻接梯度 = f_PBF_Gradient(vec_sub(x, 位置[k]), 邻接索引);
					

					uint32 邻接粒子物理ID = cCP.m_物理体ID[k];
					uint8 因子 = !(uint8(邻接粒子物理ID == 当前粒子物理ID) & 当前粒子物理类型);
					vec_multi(&邻接梯度, (-质量[当前ID] / 密度0) * 因子);

					//printf(" 位移步进 邻接索引.m_粒子光滑半径 =  %f \n", 邻接索引.m_粒子光滑半径);
					//printf(" 位移步进 -质量[当前ID] =  %f \n", -质量[当前ID]);
					//printf(" 位移步进 密度0 =  %f \n", 密度0);
					//printf(" 位移步进 x =  %f \n", 邻接梯度.x);
					//printf(" 位移步进 y =  %f \n", 邻接梯度.y);
					//printf(" 位移步进 z =  %f \n", 邻接梯度.z);
					vec_multi(&邻接梯度, (拉格朗日[当前ID] + 拉格朗日[k]));
					vec_sub(&位移步进, 邻接梯度);
					//printf("邻接梯度 =  %f %f %f \n", 邻接梯度.x, 邻接梯度.y, 邻接梯度.z);
					//printf("位移步进 =  %f %f %f \n", 位移步进.x, 位移步进.y, 位移步进.z);
				}
			}
		}
	}

	//vec_add(&位置量变化[当前ID], 位移步进);
	位置量变化[当前ID] = 位移步进;
	
}

__global__ void f_PositionUpdate(
	vec3* 位置,
	vec3* 位置量变化,
	S_DomProp			域属性
	//float32* 密度,
	//float32* 质量,
	//float32* 拉格朗日
	//S_cuAdjacencyIndex* 邻接索引
	
) {
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;

	//f_SolvingDensity()
	vec_add(&位置[当前ID], 位置量变化[当前ID]);

	位置[当前ID] = vec_Max(位置[当前ID], vec_add(域属性.m_域世界偏移, 域属性.m_单元大小*2));
	位置[当前ID] = vec_Min(位置[当前ID], vec_add(域属性.m_域世界偏移, vec_sub(域属性.m_世界大小, 域属性.m_单元大小*2)));

	位置量变化[当前ID] = { 0,0,0 };
}

__global__ void f_updateSpeed(
	vec3* 位置,
	vec3* 上一帧坐标,
	vec3* 速度,
	vec3* 加速度,
	float32* 质量
	//float32* 密度,
	//float32 粒子光滑半径,
	//S_cuAdjacencyIndex	邻接索引,
	//S_DomProp			域属性
) {
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;
	vec3 x = 位置[当前ID];

	//上一帧坐标[id] = 速度[id];
	if (质量[当前ID] != 0.0f) {
		vec3 s = vec_sub(x, 上一帧坐标[当前ID]);
		//printf("s速度s =  %f %f %f \n", s.x, s.y, s.z);
		//速度[当前ID] = vec_multi(s, 1.0 / 0.0015);
		//s = vec_multi(s, 1.0 / DEF_帧率);
		s = vec_multi(s, DEF_帧率);
		速度[当前ID] = s;
		加速度[当前ID] = s;
		//printf("速度 =  %f %f %f \n", 速度[当前ID].x, 速度[当前ID].y, 速度[当前ID].z);
		//printf(" key y =  %f", v.y);
		//printf(" key z =  %f", v.z);
		//速度[id] = 速度[id] + 加速度[id] * 更新步长;
		//位置[id] = 位置[id] + 速度[id] * 更新步长;
		//位置[id] = vec_add(位置[id], {0.0f,-0.01f,0.0f});
	}
}


__global__ void f_updateViscosity(
	vec3* 位置,
	//vec3* 上一帧坐标,
	vec3* 速度,
	vec3* 加速度,
	float32* 质量,
	float32* 密度,

	S_cuAdjacencyIndex	邻接索引,
	S_DomProp			域属性
) {
	uint64 当前ID = blockIdx.x * blockDim.x + threadIdx.x;
	vec3 x = 位置[当前ID];
	vec3 s = 加速度[当前ID];

	uvec3 coord = f_浮点坐标转索引坐标(x, 域属性.m_域世界偏移, 域属性.m_单元大小);
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 搜索坐标;
				搜索坐标.x = coord.x + ox;
				搜索坐标.y = coord.y + oy;
				搜索坐标.z = coord.z + oz;

				S_DomUnitIndex 单元索引 = f_getUnitParticleIndex(邻接索引, 域属性, 搜索坐标);
				//printf(" f_ProjectionDelta 单元索引.num =  %i \n", 单元索引.num);
				for (uint8 j = 0; j < 单元索引.num; ++j) {
					uint64 k = 单元索引.index[j];
					
					float32 w = 邻接索引.m_Viscosity * (质量[k] / (密度[k]));
					//printf("邻接索引.m_Viscosity =  %f %f %f %f \n", 邻接索引.m_Viscosity, 质量[k], 密度[k], w);
					//printf("速度[k] =  %f %f %f %llu \n", 速度[k].x, 速度[k].y, 速度[k].z, k);
					//printf("加速度[k] =  %f %f %f %llu \n", 加速度[k].x, 加速度[k].y, 加速度[k].z, k);
					vec3 v = vec_sub(s, 加速度[k]);
					//printf("v =  %f %f %f %f \n", v.x, v.y, v.z, w);
					vec_multi(&v, w);
					w = f_PBF_DensityWeights(vec_sub(x, 位置[k]), 邻接索引.m_K, 邻接索引.m_粒子光滑半径);
					//printf("v =  %f %f %f %f \n", v.x, v.y, v.z, w);
					vec_multi(&v, w);
					
					vec_sub(&速度[当前ID], v);
					//printf("速度[当前ID] =  %f %f %f \n", 速度[当前ID].x, 速度[当前ID].y, 速度[当前ID].z);
				}
			}
		}
	}
}



void f_CUDA_计算施加外力和速度叠加(
	S_PBF_Gpu& pdf, 
	vec3 重力, 
	vec2 步长范围, 
	float32& 更新步长, 
	uint64 粒子数, 
	S_DomProp& 域属性, 
	bool 是否拷回CPU)
{
	dim3 gridDim = { 1,1,1 };
	dim3 blockDim = { 1,1,1 };

	if (粒子数 > pdf.m_线程块大小) {
		blockDim.x = pdf.m_线程块大小;
		gridDim.x = 粒子数 / pdf.m_线程块大小;
	}
	else {
		blockDim.x = 粒子数;
	}


	f_forceand<<<gridDim, blockDim>>>(pdf.m_cu位置, pdf.m_cu前一帧位置, pdf.m_cu速度, pdf.m_cu加速度, pdf.m_cu外力, pdf.m_cu质量, 重力);
	cudaDeviceSynchronize();
}

void f_GPU_添加外力(S_PBF_Gpu& pbf, S_PBF粒子发射& body) {
	uint32 num = body.m_体素->m_坐标.size();

	if (body.m_偏移 > body.m_Count) {
		body.m_偏移 = body.m_Offset;
	}
	
	if (body.m_偏移 + num > body.m_Count) {
		num = body.m_Count - body.m_偏移;
		body.m_偏移 = 0;
	}

	//cudaMemcpy(&(pbf.m_cu移动变化量[body.m_偏移]), body.m_体素->m_法线.data(), num * sizeof(vec3), cudaMemcpyKind::cudaMemcpyHostToDevice);
	//body.m_体素->m_法线 这里要把VEC4数据改成VEC3
	cudaMemcpy(&(pbf.m_cu外力[body.m_偏移]), body.m_体素->m_法线.data(), num * sizeof(vec3), cudaMemcpyKind::cudaMemcpyHostToDevice);
	//body.m_偏移 += num;
}

void f_GPU_置粒子位置(S_PBF_Gpu& pbf, S_PBF粒子发射& body) {
	uint32 num = body.m_体素->m_坐标.size();

	if (body.m_偏移 > body.m_Count) {
		body.m_偏移 = body.m_Offset;
	}

	if (body.m_偏移 + num > body.m_Count) {
		num = body.m_Count - body.m_偏移;
		body.m_偏移 = 0;
	}
	
	cudaMemcpy(&(pbf.m_cu位置[body.m_偏移]), body.m_体素->m_坐标.data(), num * sizeof(vec3), cudaMemcpyKind::cudaMemcpyHostToDevice);

	body.m_偏移 += num;
}

void f_GPU_坐标钳到域内(S_PBF_Gpu& pbf, S_DomProp& 域属性, uint64 粒子数) {
	dim3 gridDim = { 1,1,1 };
	dim3 blockDim = { 1,1,1 };

	if (粒子数 > pbf.m_线程块大小) {
		blockDim.x = pbf.m_线程块大小;
		gridDim.x = 粒子数 / pbf.m_线程块大小;
	}
	else {
		blockDim.x = 粒子数;
	}

	f_Clampdown << <gridDim, blockDim >> > (pbf.m_cu位置, 域属性);
}

void f_CUDA_计算投影约束(S_PBF_Gpu& pbf, S_cuAdjacencyIndex& 邻接索引, S_DomProp& 域属性, uint64 粒子数) {
	dim3 gridDim = { 1,1,1 };
	dim3 blockDim = { 1,1,1 };

	if (粒子数 > pbf.m_线程块大小) {
		blockDim.x = pbf.m_线程块大小;
		gridDim.x = 粒子数 / pbf.m_线程块大小;
	}
	else {
		blockDim.x = 粒子数;
	}
	//m_权重0 = pbf.m_Weight0;

	S_CP_parameter cp = { pbf.m_cu位置 , pbf.m_cu物理体ID, pbf.m_cu物理体类型 };
	//std::cout << "f_SolvingDensity" << std::endl;
	f_SolvingDensity <<<gridDim, blockDim >>> (cp, pbf.m_cu质量, pbf.m_cu密度, pbf.m_cu拉格朗日, 1000.0, pbf.m_Weight0, 邻接索引, 域属性);
	//std::cout << "f_ProjectionDelta" << std::endl;
	f_ProjectionDelta <<<gridDim, blockDim >>> (cp, pbf.m_cu移动变化量, pbf.m_cu质量, pbf.m_cu拉格朗日, 1000.0, 邻接索引, 域属性);
	//std::cout << "f_PositionUpdate" << std::endl;
	f_PositionUpdate <<<gridDim, blockDim >>> (pbf.m_cu位置, pbf.m_cu移动变化量, 域属性);
	//std::cout << "f_updateSpeed" << std::endl;
	//f_updateSpeed << <gridDim, blockDim >> > (pbf.m_cu位置, pbf.m_cu前一帧位置, pbf.m_cu速度, pbf.m_cu加速度, pbf.m_cu质量, 0.001);
	//std::cout << "end" << std::endl;
	//cudaDeviceSynchronize();
}

void f_CUDA_更新粒子速度(S_PBF_Gpu& pbf, S_cuAdjacencyIndex& 邻接索引, S_DomProp& 域属性, uint64 粒子数) {
	dim3 gridDim = { 1,1,1 };
	dim3 blockDim = { 1,1,1 };

	if (粒子数 > pbf.m_线程块大小) {
		blockDim.x = pbf.m_线程块大小;
		gridDim.x = 粒子数 / pbf.m_线程块大小;
	}
	else {
		blockDim.x = 粒子数;
	}
	//f_PositionUpdate << <gridDim, blockDim >> > (pbf.m_cu位置, pbf.m_cu移动变化量, 域属性);
	//std::cout << "f_updateSpeed" << std::endl;
	f_updateSpeed << <gridDim, blockDim >> > (pbf.m_cu位置, pbf.m_cu前一帧位置, pbf.m_cu速度, pbf.m_cu加速度, pbf.m_cu质量);
	//f_updateSpeed << <gridDim, blockDim >> > (pbf.m_cu位置, pbf.m_cu前一帧位置, pbf.m_cu速度, pbf.m_cu加速度, pbf.m_cu质量, pbf.m_cu密度, 邻接索引, 域属性);
	f_updateViscosity <<<gridDim, blockDim >>> (pbf.m_cu位置, pbf.m_cu速度, pbf.m_cu加速度, pbf.m_cu质量, pbf.m_cu密度, 邻接索引, 域属性);
	//std::cout << "end" << std::endl;
	cudaDeviceSynchronize();
}


void f_CUDA_get坐标范围粒子(S_PBF_Gpu& GPU数据, S_cuAdjacencyIndex& 邻接索引, S_DomProp& dom, vec3 x, std::vector<vec3>& particle) {
	vec3* cuP;
	uint32* cuNum;

	COUDA_CALL(cudaMalloc((void**)&(cuP), sizeof(vec3) * 64));
	COUDA_CALL(cudaMalloc((void**)&(cuNum), sizeof(uint32)));


	//x = vec_Max(x, { 0,0,0 });
	//x = vec_Min(x, dom.m_世界大小 * 0.5f);

	uvec3 域坐标 = f_浮点坐标转索引坐标(x, dom.m_域世界偏移, dom.m_单元大小);
	dim3 gridDim = { 1,1,1 };
	dim3 blockDim = { 1,1,1 };



	uint32 offset = 0;
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 c;
				c.x = 域坐标.x + ox;
				c.y = 域坐标.y + oy;
				c.z = 域坐标.z + oz;

				uint32 num = 0;
				std::cout << "x:" << c.x << " y:" << c.y << " z:" << std::endl;
				f_getUnitParticle << <gridDim, blockDim >> > (邻接索引, dom, c, GPU数据.m_cu位置, cuP, cuNum);
				cudaMemcpy(&num, cuNum, sizeof(uint32), cudaMemcpyKind::cudaMemcpyDeviceToHost);

				if (num) {
					particle.resize(num + offset);
					cudaMemcpy(&(particle.data()[offset]), cuP, num * sizeof(vec3), cudaMemcpyKind::cudaMemcpyDeviceToHost);
					//vec3* p = particle.data();
					offset += num;
				}
			}
		}
	}

	cudaFree(cuP);
	cudaFree(cuNum);

	cuP = nullptr;
	cuNum = nullptr;
}









